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Abstract 

Robust heteroclinic networks are invariant sets that can appear as attractors in 
symmetricany coupled or otherwise constrained dynamical systems. These networks 
may have a very complicated structure that is poorly understood and determined to 
a large extent by the constraints and dimension of the system. As these networks are 
of great interest as dynamical models of biological and cognitive processes, it is useful 
to understand how particular graphs can be realised as robust heteroclinic networks 
that are attracting. This paper presents two methods of realizing arbitrarily complex 
directed graphs as robust heteroclinic networks for flows generated by ODEs — we say 
the ODEs realise the graphs as heteroclinic networks between equilibria that represent 
the vertices. Suppose we have a directed graph on Uy vertices with rig edges. The 
"simplex realisation" embeds the graph as an invariant set of a flow on an (n^, — 1)- 
simplex. This method realises the graph as long as it is one- and two-cycle free. The 
"cylinder realisation" embeds a graph as an invariant set of a flow on a (rie + 1)- 
dimensional space. This method realises the graph as long as it is one-cycle free. In 
both cases we find the graph as an invariant set within an attractor, and discuss some 
illustrative examples, including the influence of noise and parameters on the dynamics. 
In particular we show that the resulting heteroclinic network may or may not display 
"memory" of the vertices visited. 

1 Introduction 



Heteroclinic network attractors are an interesting nontrivial form of invariant set for nonlin- 
ear dynamics where a typical trajectory recurrently approaches a number of different unstable 
(saddle) equilibrium states. They have been found to appear robustly in various applications 
[27, 28, 20] and have been investigated by several groups as an approach to model and un- 
derstand a number of types of collective neural dynamics. In particular, the networks seem 
to be able to model a diverse range of systems, from sequence generation (deterministic or 
random) [36] and finite-state computation [39, 5], to aspects of neural function [35, 14] such 
as binocular rivalry [7]. Close dynamical relatives are found in models of phase- resetting 
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oscillators where "unstable attractor networks" are a limiting case of heteroclinic networks 
that can appear in non-invertiblc scmiflows [11, 15, 30]. 

Until now, the main problem with using such an approach to model, for example, a 
specific switching process, has been that coupled systems can possess quite complicated 
robust heteroclinic networks. It is highly nontrivial to determine whether a coupled cell has 
a heteroclinic network attractor, and even when one can find them, they may have very high 
complexity on increasing the number of coupled cells. For example, heteroclinic networks 
the structure of "odd graphs" [8] can be found in systems of = 2A; + 1 oscillators, and 
heteroclinic ratchets can be found in less symmetric coupled systems [22]. Up to now there 
seems to be no way to easily design dynamical systems and/or coupling between them that 
will realise a given directed network as a robust heteroclinic attractor. The closest approaches 
we are aware of are [33] and [1] where heteroclinic networks are designed in small numbers 
of coupled cells by a process of "inflation" of cells within a smaller network. However, the 
networks must be reducible to combinations of simpler cyclic structures. 

The main contribution of this paper is to give two explicit methods of how to design 
a coupled cell system that realises any given directed graph as a heteroclinic network. In 
doing so, we suggest possible robust ways to embed a finite-state discrete state computational 
system into a system of autonomous coupled dynamical systems. The method in Section 2.1 
(which we call the simplex realisation) realises any graph (that is one- and two-cycle free) 
as a heteroclinic network for a cubic polynomial vector field on a simplex. The method in 
Section 2.2 (which we call the cyhnder reahsation) reahses any graph that is one-cycle free 
based around a vector field on a line coupled to a number of "transition modes" . We do this 
by constructing a "coupled cell system" that in both cases is "all-to-all" coupled of a small 
number of cell types, but where the coupling is quite inhomogeneous. Figure 1 illustrates 
the realisation methods. 

We analyse both of these realization methods and illustrate their application to some 
example graphs in Section 3. In all but the simplest cases the embedded network will be 
part of a larger chain-recurrent set that includes additional "induced vertices" and separa- 
trices between different connection sets; nonetheless we find in numerical simulations that 
these induced vertices do not unduly infiuence the the behaviour of the network, although 
we have no proof of this. In Section 4 we consider long-term statistics of the behaviour near 
the network under the infiuence of noise fiuctuations. We see that the transition behaviour 
between vertices is in certain cases well-modelled as a Markov process with transition de- 
pending only on the current vertex. However, the presence of "lift-off" [3] of trajectories 
can lead to short-term memory — i.e. cases where the transition probability depends not only 
on the current vertex but also on the past few vertices visited; we lose the Markov prop- 
erty. This is observed in numerical simulations and analysed discussed in Section 5. Finally, 
Section 6 finishes with a discussion of issues around the proposed methods of realization of 
heteroclinic cycles, and we give general thoughts on how these constructions may be related 
to applications. 
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Simplex realisation Directed graph Cylinder realisation 



Figure 1: The two realisation methods described in this paper for a coupled cell network 
to realise a specific directed graph with vertices and rie edges as a heteroclinic network. 
The "simplex" realisation described in Section 2.1 uses identical cells Xi with an inho- 
mogeneous coupling while the "cylinder" realisation described in Section 2.2 has a special 
cell (p) that remembers the last vertex visited while the other Uf. cells (t/j) are identical and 
have inhomogeneous coupling to provide dynamical connections between the states. The 
dashed lines indicate that the nature of coupling between the Ui is of a different nature to 
that between that Ui and p. The names are due to invariant structures in phase space rather 
than their coupling. 

2 How to design heteroclinic networks with a given 
directed graph structure 

Suppose that Q = (y,S) is a directed graph between a finite set of vertices V = {o,}"^]^ 
with directed edges £ = {ti}^^^. Let a(Cj) denote the starting vertex and oj{ti) the finishing 
vertex of an edge Cj. We will write Ajk to be the adjacency matrix for the graph, i.e. 

^ _ f 1 if there is an i with D-,- = a{te), D/c = uj{te) 
■'^ y otherwise. 

Now consider the ODE on M"'' describing a set of coupled cells 

dx ■ 

= fi{xi,...,Xn)+CWi{t) (1) 

where Xi G M"' for i = l,...,n, Wi represents additive i.i.d. noise and C ^ the noise 
amplitude. We assume that / = is sufficiently smooth (at least so that 

linearizations vary robustly). We say an invariant set X C M""' for the fiow induced by (1) 
realises the graph ^ as a heteroclinic network if X is an embedding of Q such that each 
vertex Oj is mapped onto an equilibrium for the fiow (contained in X), and each edge 
is mapped onto (possibly one of a set of) heteroclinic connections from a{ti) to a;(ej), also 
contained in X. 
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It is well known that a realisation may be robust if the form of / is constrained — see for 
example [27, 17, 6]. The network may also be attracting if certain eigenvalue conditions are 
satisfied. Note that (a) robustness implies necessarily that all equilibria will be hyperbolic 
and (b) both of our constructions give robust realisations in that there is a a symmetry 
group generated by order two elements, and the constructions are robust to perturbations 
that respect this symmetry group. 

In the remainder of this section, we present two methods for creating heteroclinic networks 
in a coupled cell network with d = 1 dimensional dynamics within each cell; the coupling 
structures are sketched in Figure 1. The first is a construction where the vertices of the 
graph are embedded as vertices of an (n^, — l)-simplex (note that an n-simplex has n + 1 
vertices) — we call this the "simplex reahsation" . Many examples of heterochnic cycles and 
networks in the literature are of this type, especially the "winnerless competition/Lotka- 
Volterra" type models. See, for examples, [18, 23, 33]. The second method, inspired by 
[7, Figure 11] has all vertices embedded on an invariant line on the centreline of an rig + 1 
dimensional cyhnder — we call this the "cyhnder reahsation" . We will be interested in the 
long time (asymptotic) behaviour of the system in the presence of noise and inhomogeneities, 
where the heteroclinic network will not be an exact solution but it will continue to give a 
very clear framework in which one can discuss the dynamics of nearby solutions. 



2.1 Simplex realisation 

The first construction proceeds as follows. Suppose that ^ is a directed graph with vertices 
and consider the stochastically forced vector field on a; e M.""" for the simplex realisation as 
follows: 

^ - Xj (^1 - \x\' + J2 j + C^^^(^) (2) 

with Xj G M for j = 1, . . . , n„ and |a;p = x^. The wj represent i.i.d. white noise modulated 
by < C <^ 1. Note that in the absence of noise, the system has Zg" symmetry, where each 
of the Z2 subgroups is generated by a reflection in the Xjth coordinate. This symmetry 
implies that each coordinate hyperplane is invariant under the flow (for compactness, we use 
the convention that a set is said to be invariant if it is invariant for the noise- free system). 

We say a graph Q is n-cycle free if is contains no directed loops of length n. In the 
following proposition we will consider graphs where there are no edges with a{ti) — uj{ti) 
(no one-cycles) and no pairs of edges e^j with a{ti) — uj{tj) and a{tj) — uj{ti) (no two- 
cycles). Note that Q is n-cycle free if and only if all diagonal entries of A" are zero, where 
A is the adjacency matrix. 

Proposition 1 Suppose that G is one- cycle and two-cycle free. Then Uij can be chosen so 

that the system (2) for C = realises the graph Q as a heteroclinic network. Moreover, this 
realisation is robust to perturbations that respect symmetry given by reflection in the 
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coordinate planes} 

Proof: The system (2) has equihbria on each of the fiy coordinate axes, 

a = (o,...,o,i,o,...,o) 

where the 1 is in the kth position and k = 1, • • • , n^. We associate with the vertex 0^ of 
the graph Q. The eigenvalues of are —2 + Sukk in the radial {xk) direction, and akj in the 
Xj direction {j ^ k). We choose the coefficients aij so that ajj — 0, and 

Gij > if Aij = 1 while a^j < if Aji = 0, {i ^ j). 

If the graph Q contains an edge (e;) with a{ti) — and uj{ti) = (and hence does not 
contain an edge from Oj to Oj), then Aij = 1, hence aij > and so the equilibrium ^i will 
be unstable in the Xj direction. Similarly, Aji = 0, thus aji < and so C,j will be stable in 
the Xi direction. Restricted to the two-dimensional subspace {xi,Xj), the resulting flow and 
connecting heteroclinic orbit is shown in Figure 2(a). 

Because the network is one- and two-cycle free we have Aij Aji — for all If this 
condition were broken, for instance if Aij = Aji = 1, then $^i would be unstable to perturba- 
tions in the Xj direction, and ^j would be unstable to perturbations in the Xj direction. This 
would result in a flow as shown in Figure 2(b), and a new (stable) equilibrium is created in 
the {xi,Xj) plane. QED 



This system also contains a number of other equilibria; for example, the equilibrium at 
the origin, which is unstable (all eigenvalues are equal to 1). The coefficients aij that are 
non-zero can be thought of as weights on the graph; the magnitude of a^- > affects the 
expanding eigenvalue at while Uji < affects the contracting eigenvalue at ^j. In later 
sections we will parameterize the dependence by making the a function of a matrix of 
weights Wij > 0. 

This construction gives a realisation of the graph Q that may be embedded in a larger 
heteroclinic network that could realise a larger graph; it is also possible that the network is 
not attracting. In general, computation of the stability properties of a heteroclinic network 
is quite involved, see for example [24, 25]. We conjecture that a sufficient (but by no means 
necessary) condition for stability of the heteroclinic network is that all the contracting eigen- 
values are greater than all the expanding eigenvalues. This can be achieved in this example 
by choosing the aij appropriately. The computation of detailed stability conditions for the 
network is beyond the scope of this paper, but we find in our numerical examples that the 
networks can be easily observed and hence appear to be attracting. 

^In particular, the robustness to perturbations implies that the construction is robust to choice of param- 
eter values. 



5 




Figure 2: (a) Schematic diagram showing a connecting orbit between equihbria and in 
the (xj, Xj)-plane for the simplex reahsation. (b) The lack of connection in this invariant 
subspace if aij > and aji > 0; there is an additional equilibrium that is a sink within the 
subspace. Equilibria are shown with dots; the asymptotic dynamics for the full system is on 
topological {uy — l)-simplex within an attracting invariant sphere. 

2.2 Cylinder network 

The second construction also proceeds directly from the graph structure. Suppose that Q is 
a directed graph with with He edges; let us define the stochastically forced vector field for 
the cylinder realisation on {y,p) G M""^"'"-'^ by 

^ = -yjG,{y,p) + Cw,{t) (3) 

dp 

— = - sin(27rp) + Fj{y,p) + Cwo{t) 

for j = 1, . . . where Fj and Gj are smooth functions and are even in each of the yj, 
and Fj{0,p) = for any p. As before, the wj are i.i.d. white noise processes with amplitude 
< C < 1. 

This system (with C = 0) has an invariant line yi = y2 = ' ' ' = Vn^ = parametrized by 
p, on which there are equilibria at p = n, n G Z, yj = 0. We denote the equilibrium 

efc = (0,...,0,fc) 

for k = 1, . . .Uy. The invariant line is contained in each of the invariant planes i = 1, ■ ■ ■ , ng: 

Pr.= {{y,p) : yj = o if j^i}. 

We say the variable is activated (i.e. it is non-zero and all other yj remain zero) in P^, 
on the connection from a{te) to u;(e^). The following result holds for any graph as long as 
there are no edges Cj with a{ci) = uj{ti) — it constructs a network that lies in a cylindrical 
neighbourhood of the invariant line — hence the name of the network. 
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Proposition 2 Suppose that Q is one-cycle free. Then there are smooth functions Fj andGj 

such that the system (3) for C = realises the graph Q as a heteroclinic network. Moreover, 
this realisation is robust to perturbations that respect IJ^" symmetry given by yj ^ —yj. 

Proof: We consider the following specific choices for Fj and Gj: 
Fj{y,p) = ^y|tanh(a;(ej) -p), 



Gj{y,p) = 



(4) 



faj{p) = LajSech^ {kaj{p - ai^j))) , 

fwj{p) = L^jSech^ {koj j{p - ^(fj))) ■ 

The quantities Laj, L^^j, kaj, ki^j, Ki > are parameters that will be chosen appropriately so 
that there is a robust heteroclinic connection in corresponding to the edge e^. Note that 
Ki is a "mutual inhibition" parameter between the various yj. 

If Lae = L^i = then, within the plane Pi, the system has y^-nuUclines at y^ = 0, 
y^ = 1/2 and yi = 5/2. The additional terms when L^^, L^^ ^ make a{t() unstable in the 
direction near p — Q;(e^). Constraints on the parameters can be understood by considering 
the geometry of the heteroclinic connection from a{tc) to uj{te), which lies in P^. 

• Firstly, we require that ye < for all ye at p — uj{te), which is achieved if 

y^-l) -l-/a^Me,)) + L,, >0, Vy 

which is satisfied if Li^e > 1 and faiij^i^)) is sufficiently small. 

• Secondly, we require that y^ > when p — a{t() and < ye < yu for some y„ > 1/2. 
This is achieved if 

2/^ - ^ j -l-Lae + fc.e{a{te)) < 0, < y < y„ 

which is satisfied if Lae > 9/16 and fue{<^{^e)) is sufficiently small. 

Each of the equihbria has eigenvalues as follows: firstly, all equilibria have a 'radial' 
eigenvalue in the p direction which we label = — 27r. Now, suppose there is an edge Cj 
with a{e.i) — 0^, then the eigenvalue at in the y^ direction is 

9 

Id 

Similarly, if there is an edge with uj{tj) — t)jt, then the eigenvalue at in the yj direction 
is 

9 

C-kj ~ ~Yg faj{k) ~ Li^j. 



If the edge neither starts nor ends at 0^ (i.e. a{ti) ^ 0^ and a;(e^) ^ 0^) then the eigenvalue 
at in the yi direction will be 

9 

Sufficient conditions for the existence of the desired connections are that eki > 0, for all i for 
which a{Ci) ~ fk (which gives the condition L^i > 9/16, as before) and Ckj < 0, for all j for 
which uj{Cj) = Dfc. 

Hence we can choose parameters Laj, L^^j, kaj and k^^j so that faj{p) and fujj{p) are small 
enough when p is far from a{zj) and u!{tj) respectively. Alternatively, we can set the values 
of the contracting and expanding eigenvalues (subject to some constraints; see below), and 
then choose parameters L^j, Li^j, kaj and ki^j so that the network has this set of eigenvalues. 
QED 

Figure 3 shows the structure of the connections in phase space using this construction. 
One implication of the requirement that L^^j > 1 is that for these parameters, the contracting 
eigenvalues will always be less than —25/16. The final term in the yj equations, Ki ^i^j yf, 
is a "mutual inhibition" term that is included so that additional stable equilibria away from 
the coordinate planes are suppressed. We will typically choose = 1. Similarly to the 
simplex construction, the realisation of Q may be part of a larger heteroclinic network and 
the network may not be an attractor unless care is taken to ensure that the expanding 
eigenvalues are weak enough. A necessary (but not sufficient) condition for asymptotic 
stability of the heteroclinic network is that tk£ < 0. 

Within the two-dimensional subspace Pg parameterized by (p, yi) we show nuUclines of 
the system with typical parameter values in Figure 4, in the case of a{te) — 1 and uj{te) — 4. 
The blue lines show the p nuUchnes while the red lines show the yg, nuUchnes. Note that 
within the subspace yj = 0, the equilibria are sinks. Also, note that there are additional 
equilibria at p = A; + |, which are unstable in this one-dimensional subspace. Trajectories 
which start near move towards ^4. 

3 Examples of networks realised as heteroclinic at trac- 
tors 

We now consider two illustrative examples; the first (the decision graph) uses the simplex 
realisation and Proposition 1 while the second (the Petersen graph) uses the cylinder realisa- 
tion and Proposition 2. The decision graph could be realised using the cylinder realisation, 
but the Petersen graph cannot be realised using the simplex realisation as it contains cycles 
of order two. 
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Figure 3: Schematic diagram showing the structure of part of a cyhnder reahsation: the 
vertices hue up along the p-axis while the {yj,p) planes contain connections between vertices. 
In this case there are connections ^3 — )■ ^4 in Pi, ^1 — )■ ^3 in P2 and ^1 — )■ ^2 in Ps- Note that 
the planes Pi are all orthogonal. 
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Figure 4: The figure shows the nullclines and vector field for a typical two-dimensional 
subspace of the cylinder realisation; with a connection from p = 1 to p = 4. Arrows in the 
vector field are all scaled to have the same length for clarity. In this case we have a connection 
from the saddle ,^1 at p = 1 to the sink ^4 at p = 4. Note that there are additional sinks at 
p = 2, 3, 5 etc that are surrounded by bounded basins of attraction. 
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Figure 5: The "decision network" graph for the adjacency matrix corresponding to the 
weights in (5). 



3.1 Decision graph 

Consider the "decision tree with reset" shown in Figure 5, which for convenience we call the 
"decision graph" . One can think of this graph ^ as a two-level binary decision tree followed 
by a reset back to the first (root) vertex of the tree. This graph has = 8 vertices, = 11 
edges and an adjacency matrix corresponding to the non-zero elements of 
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If we set parameters 



and define 



0.99, a' = 0.98, a = 0.2, fi = 0.6, C = 10" 



(6) 



-a 



flWij 







otherwise 



then it is simple to check that the requirements of Proposition 1 are satisfied and (2) will 
realise a heteroclinic network with the structure of the graph shown in Figure 5. 

One can visualise the dynamics on the realisation of this graph by projecting the trajec- 
tory onto two observables that order the vertices into a ring. More precisely, for a trajectory 
x{t) we define the complex observable 



R{t) 



n 

-E 



Xj. exp 



tTC- 



2{k-l] 



k=l 



(7) 
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so that the vertices of the graph are projected onto the vertices of a regular n^-gon on 
the unit circle. Figure 6 shows a projection and timeseries for the dynamics of the system 
corresponding to the graph in Figure 5. 

Proposition 1 means that the graph structure is unaffected by the exact values of the 
parameters used. However, the statistics of the residence times and the transition proba- 
bihties at decision points (including the possibihty of memory) are strongly affected by the 
parameter values. We explore some of these effects in Section 4 and apply our results to 
investigate this graph in Section 5. 

3.2 Petersen graph 

As an example with a somewhat different structure (and one that is nontrivially but highly 
connected), wc now consider a realisation of the Petersen graph (shown in Figure 7) using 
the cylinder realisation. Each edge can be thought of as a pair of directed edges, one in each 
direction. 

This graph has been studied extensively as one of the simplest examples of a graph 
with certain nontrivial colouring properties; it has also been found to organize heteroclinic 
networks in five globally coupled phase oscillators and in systems of five delay pulse-coupled 
oscillators, where it has been proposed for computational purposes [5, 31]. 

Figure 8 shows a simulation of a reahsation of this using the cyhnder reahsation (3,4) 
and adjacency matrix 

110 0" 

1 

1 1 

1 1 

1 1 , . 
110' 
1 1 
1 1 
10 110 
1 1 1 _ 

The other parameters used for the simulation in Figure 8 are 

L^j = 1.4376, L^j = 1.5625, kaj = 2.017, k^j = 0.4705, = 1, ( = 10"^ (9) 

Note that by varying kaj and k^ij with j we can vary the expanding and contracting eigen- 
values on the jth connection. For simplicity, in the example we set them to be identical. 

3.3 Further examples 

We briefly highlight two further graphs that illustrate different ways in which vertices can 
connect cycles of the same type within a network. Figure 9 illustrates two networks each 
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Figure 6: An example simulation of the decision graph (shown in Figure 5) using the simplex 
realisation (2), using the weights (5). Top: time series of xj; the red (resp. black) dots on the 
lowest frame indicate times where the trajectory enters (resp. leaves) an "epoch" of being 
close to one of the saddles. Below: projection onto the real and imaginary parts of order 
parameter for the system showing the location of the equilibria. Simulations clearly show an 
attracting heteroclinic network with the structure as in Figure 5, projected into the plane 
using the complex observable R(t) (equation (7)). Note that presence of additional saddles 
can be inferred from rare excursions away from the one-dimensional connections. 
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Figure 7: The Petersen graph; this can be viewed as a (one-cycle free) directed graph by 
considering each edge as a pair of directed edges in both directions. 




50 100 150 200 250 300 
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Figure 8: Time series for a cyhnder reahsation of the Petersen graph (Figure 7). The p- 
dynamics (top) can be observed to wander between the vertices of the Petersen graph, only 
making transitions corresponding to edges in the graph, while each of the components i/i for 
i = 1, . . . , 15 (bottom) become non-zero only during a transition along the ith edge. The 
presence of weak noise causes the dynamics to wander around the network. Note that all the 
saddles corresponding to vertices on the network have three dimensional unstable manifolds. 
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Figure 9: (a) A "Bowtie" graph of two connected cycles and (b) the "Kirk-Silber" [23] 
network that shows a competition of two length three cycles with a shared edge. Note that 
in both cases there are two cycles of length three; in (b) there is a shared edge and two 
shared vertices while while in (a) there is only a shared vertex. 

containing two nontrivial cycles of length three. In each case the network can in principle be 
realised using either "simplex" or "cylinder" realisations. Note that although the decision 
graph has alternative routes around the graph, all of these share at least one edge; similarly 
the cycles on the Kirk-Silber network share an edge, while those in the Bowtie graph do not. 

4 Statistical properties of trajectories near a realised 
network 

Propositions 1 and 2 do not give unique ways to realise a given graph ^ as a heteroclinic 
attractor - rather, they give open sets of functions that give the appropriate embedding. 
These alternative realisations can however display different statistical properties of trajecto- 
ries near the network in the presence of low noise (or more generally, time- dependent inputs) 
after transients have decayed. In this section, we adapt and generalise some of the analysis of 
Stone, Armbruster and colleagues [38, 37, 3] to understand some basic statistical properties 
of the noise-induced itineraries of trajectories near the network. 

We fix on a size h > for neighbourhoods of the equilibria and say a solution trajectory 
x{t) is close to an equilibrium (corresponding to a vertex of the realised graph) at time 
t if \x(t) — < h. We require that — > 2h for any p q so that at any time, 
the trajectory is close to at most one equilibrium. We say x(t) remains close to C,j during 
t E [s,s + T] if it satisfies 

\x(t) — < /i for t G (s, s + T), and \x(t) — = h ioi t E {s, s + T}. 

We divide a trajectory starting at some given initial condition x and given noise path Co into 
an itinerary, i.e. the sequence of epochs: 

{{ik^Sk^Tk) : fceN} (10) 
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such that x{t) remains close to ^j^. during t e [sk,Sk + Tfc]- Note that ik G 
represent indicies of the vertices (equihbria), the times of entry < are increasing and the 
durations T^. > are all positive; for example, the red and black dots in Figure 5 illustrate 
the sequences Sk and Sk + Tk respectively, for h = 0.1. 

One can formally consider the solution of the noise-perturbed system as deterministic in 
the spirit of random dynamical systems [2], i.e. as a skew product evolution 

x{t) = M^,^) ni^ 

over the evolution 6t along the noise path. In the original system u represents a particu- 
lar Brownian paths whereas in the simulation, a) can be thought of as choice of seed for 
a random number generator. Note that satisfies the cocycle property ^t+s{x,o^) = 
^t{^s{x,uj),9s{uj)). Then x(t) remains close to if a;) — < h, which clearly is 

determined by initial condition and noise path. 

If the trajectory remains close to the attractor then it will have an infinite itinerary 
except in the (very unlikely) case that it remains close to one of the saddle equilibria for 
all time. We make a stochastic stationarity assumption that the statistical properties of 
the itinerary of typical initial conditions and typical noise trajectory are stationary and 
independent of the initial condition and details of the noise trajectory. This means we 
assume that any transients associated with the initial condition will have decayed and that 
the initial condition is distributed according to a stationary invariant probability distribution 
in the space of initial conditions x for the noise-perturbed system, and the noise path u is 
"typical" . More precisely, we make an implicit assumption that there is an ergodic probability 
measure for 6t that lifts to a "natural measure" for the skew product flow. With respect to 
this distribution, we define the probability of observing a given finite sequence of vertices 
{jk : k ^ 1, . . . ,m} as 

"^(ji) • • • ) jm) = Pi'ob( {x, uj) : ie = je for i = 1, ■ ■ ■ ,171, for the trajectory starting at {x, u)). 

(12) 

The stationarity assumption implies that one can compute this from a typical trajectory 
simply in terms of the frequency of transitions 

^(j'l, ■■■,jm) = lim 7#{0 <i <k : ii+n = jn forn = 1, . . . , m}. (13) 
This can be used to define the probability tTj of an epoch being close to ^f. 

TT,- = P(j) 

and (assuming ttj > 0) we define the transition probability between vertices and by: 

= (14) 

We ask the question: are the itinerary statistics to those of a Markov chain whose non-zero 
transitions correspond to edges of the original graph Ql Although this is possible, subtle 
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effects can appear (see [12]) wfiere transition probabilities may depend not just on tfie current 
equilibrium but also on previous equilibria visited. 

We say a noise perturbed attracting heteroclinic network is memoryless if the transition 
probabilities are independent of the prior itinerary, i.e. if the itineraries are Markov of order 
one. More precisely, a transition from to is said to be memoryless if 

^(il, • • • , im) = ViJl, . . .,jm-l)7^p,q (15) 

for any m > 2 and any sequence {jk : k = 1, . . . ,m} with jm-i = P and j„i = q. The 
system is memoryless if all possible transitions are memoryless. We expect such a system 
to be truly memoryless only in the singular limit of very low noise C ^ 0; nonetheless, later 
sections suggest that a system can be very close to memoryless in that corrections to (15) 
may be asymptotically small in (. Note also that a particular transition will necessarily be 
memoryless if TTp^q — 1, though this does not necessarily imply that all other transitions are 
memoryless. 

The transition probabilities will be affected (to a greater or lesser extent) by the eigenval- 
ues at the equilibria and by the noise level. In the following section, we outline the influence 
of these parameters. Finally, using the simplex realisation of the decision graph. Section 5 
gives some numerical examples. 

4.1 Analysis of dynamics near heteroclinic networks: Poincare 
sections and maps 

Analysis of the dynamics of trajectories near heteroclinic cycles and networks is often done 
via construction of Poincare maps to approximate the flow (see, e.g. [23]). The flow is divided 
into two parts: a local flow defined in a neighbourhood of the equilibria and a global flow 
away from the equilibria, usually constructed by linearising the flow around the heteroclinic 
orbits. The local and global flow define maps between Poincare sections usually defined close 
to equilibria; their composition gives a map which well approximates the flow close to the 
heteroclinic network. 

The detailed construction of such maps for the networks (even in the noise-free case) is 
beyond the scope of this paper. However, in later sections we apply previous results regarding 
the dynamics of trajectories near heteroclinic networks in the presence of noise, so we find it 
useful to define the appropriate Poincare sections here. Poincare sections are often defined to 
be surfaces a distance h from an equilibrium, where h is some small constant, meaning they 
are spheres. However, in order to apply results of Stone, Armbruster and others [38, 37, 3] 
we define Poincare sections as unions of codimension one surfaces as follows. We give explicit 
definitions of Poincare sections for the simplex realisation; for the cylinder realisation they 
can be defined similarly. 

Consider an equilibrium in the simplex realisation, with an incoming heteroclinic 
connections from equilibria and an outgoing heteroclinic connection towards equilibria 
Then we define (for a fixed small h): 

H'I''^ = {x e E"" : Xj- = /i, |1 - xjtl < h, \xi\<h,i^ j, k} (16) 
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Figure 10: Schematic diagram showing a Poincare sections H^^'^, ifg"*'^ cind ifg*^*'^ near the 
equihbrium ^2- Heterochnic orbits are shown with bold hnes. Large dots indicate equihbria, 
small dots indicate the intersection of the heteroclinic orbits with the Poincare sections. The 
dotted lines on if™ indicate the dividing line between regions of trajectories which will next 
travel to ,^3 or ^4, in the case that 623 > 624. 



= {xe M"" ■.xi = h,\l-Xk\<h, \xi\ <h,iy^ l,k}. (17) 

If a vertex in the graph has an incoming or outgoing degree greater than one, then there 
will be multiple incoming or outgoing Poincare sections, and so we also define 

j I 

where the unions are taken over all incoming and outgoing directions respectively. Note that 
if the trajectory spends time near it must have passed through if™ and will pass through 
if^"* at some time thereafter. Figure 10 shows a schematic diagram of the Poincare sections 
near an equilibrium ^2, with an incoming connection from ^1 and outgoing connections to ^3 
and ^4. 

In the following sections, we use the following notation for eigenvalues of equilibria. If 
has a contracting direction in the xj direction, we label the corresponding eigenvalue —Ckj- 
Similarly, if If has an expanding direction in the xi direction, we label the corresponding 
eigenvalue Cki- 

4.2 Transition probabilities between equilibria 

If has a single expanding direction towards then clearly tt^- ^ = 1 and vr^y = for k ^ I'm 
the limit of low noise. The more interesting case is when the equilibrium has two expanding 
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directions, as in the example shown schematically in Figure 10. Here, at the equilibrium ^2, 
the trajectory makes a 'choice' as to whether to next visit ,^3 or ,^4. We give here an outline of 
how to compute the probabilities of the trajectory making each choice; that is, the transition 
probabilities. The extension of these calculations to equilibria with three or more expanding 
directions is straightforward. 

Armbruster et. al [3] compute the transition probabilities for an equilibrium with two 
expanding directions, under the assumption that the incoming distribution of coordinates 
of a trajectory is approximately Gaussian and centered at zero. For the schematic example 
in Figure 10, we show this distribution of incoming trajectories, in terms of the 2:3 and 
X4 coordinates on H™, on the left-hand side of Figure 11. The proportion of trajectories 
which we expect to visit ^4 is given by the proportion of the area of the 'noise ellipse' which 
intersects the cusp. 

The computation of this area was given in [3]. The results depend on constants which 
come from the global part of the flow and are in general unknown. However, a scaling can 
be found in the limit of low noise. It was shown that if an equilibrium has expanding 
directions Xk and Xm, with expanding eigenvalues Cjk > ejm, then in the limit of low noise, 
(i.e. as C — >^ 0) the transition probabihty from to scales as 



4.3 Lift-off and memory effects for transitions between equilibria 

For some applications, it may be desirable to include memory effects into a network de- 
sign where the transition probabilities depend on the recent history of the trajectory taken 
through the network. It is possible to create memory effects in a noisy heteroclinic network 
that have "lift-off", an effect noted by Stone, Armbruster and Kirk [37, 3] that appears for 
attracting heteroclinic (but not homoclinic) networks with additive white noise. 

Lift-off is a property of the distribution of the coordinates of a trajectory as it enters 
a neighborhood of an equilibrium such that the probability distribution of coordinates of 
trajectories in a Poincare section becomes multi-modal and is a mechanism by which the 
transition probabilities can gain memory. We demonstrate the idea of lift-off in Figure 11, 
which shows the Poincare section i?™'^ for the schematic in Figure 10, and a representation 
of the distribution of the coordinates of the incoming trajectories. If no hft-off occurs, then 
the distribution of coordinates is approximately Gaussian and centered at zero, as shown 
on the left. If lift-off occurs, then the distribution of a particular coordinate is no longer 
Gaussian and may include several peaks, as shown on the right. 

The transition probabilities can be thought of as intersections of the distribution of 
incoming trajectories with the the set of points that progress to a given equilibrium at the 
next epoch. As lift-off changes this area of intersection and depends on the previously epochs, 
lift-off can affect transition probabilities and make them depend on memory. The occurrence 
of lift-off depends on the eigenvalues of the equilibria, and hence can be controlled. Note that 
a sufficient condition for lift-off not to occur is that all contracting eigenvalues are greater 
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Figure 11: Schematic diagram showing the Poincare section H'^'^ for the connection from 
to ,^2 shown in Figure 10, projected onto the X3,X4 plane. The dotted hnes indicate the 
dividing hne between regions of trajectories which travel next to .^3 or ,^4. The shaded ellipse 
represents the distribution of coordinates of trajectories as they pass through the section 
where the darker regions are visited more often, while the black dot shows the location of 
the noise-free connection. The left panel shows a case in which the incoming distribution of 
the X3 and X4 coordinates is approximately Gaussian and centered at zero. The right panel 
shows a case with lift-off in the direction and consequently a multi-modal distribution. 



than all expanding eigenvalues. 

It was shown in [37] that if the network contains a heteroclinic connection — > ^k, then 
there is the possibility of creating lift-off in the Xj direction as the trajectory exits a neigh- 
bourhood of the equilibrium ^fc. If this lift-off can be maintained until the trajectory next 
approaches an equilibrium with an unstable manifold which connects to ^j, the probability 
of visiting will be higher than if the lift-off had not occurred. 

Stone and Armburster [37] compute conditions on the eigenvalues for lift-off to occur. 
They show that, if there exists a connection — >■ ^fe, lift-off will occur at (that is, in the 
distribution of coordinates of trajectories leaving a neighbourhood of ^k) in the Xj direction 
if 

^<1, 

where Xi is the expanding direction at That is, if this condition is satisfied, then the 
distribution of the xj coordinate as the trajectory exits a neighbourhood of C,k is not centered 
about zero, but instead is centered about a point whose distance from the connection scales 
as 

in low noise limit C ~^ 0. Since the distribution of the noise is symmetric about zero, the lift- 
off has equal probability of occurring in either the positive or negative direction. This means 
that the distribution becomes multimodal. Since the global flow only scales the distribution 
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of coordinates by an amount that is order one in (, the incoming distribution to the next 
equihbrium will scale similarly. 

It is straightforwards to extend the calculations of [37] to a sequence of heteroclinic 
connections between equilibria ^ ^ Cm- That is, we already have conditions that 

the distribution of the Xj coordinate as it exits a neighbourhood of is not centered about 
zero, and we can extend the computations to give conditions that the distribution of the Xj 
coordinate as the trajectory exits a neighbourhood of is similarly not centered about zero. 

These calculations are a simple extension of those given in [37] , and it can be shown that 

^ + ^<1, (18) 

^kl ^Im 

implies that the distribution of the Xj coordinate as the trajectory exits a neighbourhood of 
is centered about a point whose distance from the connection scales as 

in low noise limit C ~^ 0. Again, since the global flow only scales coordinates by an order 
one amount, the incoming distribution of the Xj coordinate at equilibrium will then also 
not be centered at zero. Thus, if condition (18) on the eigenvalues is satisfied, and has 
an unstable manifold (of dimension two or greater) which includes the direction towards ^j, 
then the proportion of trajectories which go towards will be higher than if lift-off had not 
occurred. 

Note that lift-off occurs only if the contracting eigenvalues are small enough, that is, 
if there is not enough contraction in the Xj direction at to 'squash' the lift-off. Further 
extension of these calculations will give conditions on lift-off being maintained over longer 
sequences of equilibria, and hence the possibility of longer term memory. However, as lift- 
off requires sufficiently small contracting eigenvalues, and stability of the network requires 
sufficiently large contracting eigenvalues, we conjecture that any memory effects present must 
appear of a sequence of epochs that is smaller than the longest cycle within the network. 

5 Example: simplex realisation of the decision graph 

In this section, we numerically investigate the properties described in Section 4 for the 
decision graph realised via the simplex realisation as discussed in Section 3.1. Note that 
whilst it is comparatively easy to check for the presence of memory in solution trajectories, 
checking for the absence of memory is difficult. That is, although we can perform statistical 
tests to show that a trajectory may be zeroth order rather than first order Markov, checking 
for the absence of all long-term memory effects is hard. Whilst we expect that such long 
term memory (i.e. over sequences of more than two equilibria) is possible, it seems unlikely 
to occur in the absence of short term memory except in special cases. A detailed study 
of long-term memory is beyond the scope of this paper and so in the following sections we 
discuss only the presence or absence of short-term memory. 
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5.1 Transition probabilities in the decision graph 

First, we consider the case where we predict that the network is memory less, that is, there 
is no hft-off, and so the incoming distributions of all coordinates at each equilibria are 

Gaussians centered at zero. 

We confirm the scaling given in Section 4.2 for the transition probabilities in the decision 
graph by performing numerical integrations at various noise levels. We use parameters Cjk = 

tjk = 2, ei2 = 2 , 624 = 636 = 648 = 1-85, 623 = 635 = 647 = 1-99, 651 = 661 = 671 = 631 = 1.98, 

and integrate for noise levels C = 10"^ 10"^ 10^^ W~^, 10"^ 10"^°, 10"^^ At each noise 
level we integrate until the number of passes the trajectory makes through is 5000. 

At each noise level, we measure the number of visits the trajectory makes to the equilibria 
^4, ^6 and as a proportion of the number of visits to ^2, and respectively. In terms 
of the definitions given in Section 4, these are approximations of 

— , — , and — 

712 TTs 7r4 

respectively, using approximations (13) and k = 5000 x 4 = 20,000 with standard errors. 
These proportions are plotted against the noise level in Figure 12 on log-log axes. Best fit 
hnear regressions for each of these lines have slopes 0.072, 0.074 and 0.073 respectively. The 
expected slope for all three is — 1 — 0.076. This means that in the low noise limit C ~^ 
we predict that almost all itineraries follow the cycle 

^ ^ 6 -> 6 -> ^5 ^ 6 ^ • • • 

corresponding to selecting the most unstable expanding direction at each equilibrium. 

Note that without knowing details about the global part of the fiow, it is not possible 
to predict the precise transition probabilities for particular noise levels and eigenvalues, the 
best we can do is get the scaling for low noise. However, if a particular transition probability 
is desired, a little experimentation and alteration of noise levels can achieve this. 

5.2 Numerical examples of memory effects for the decision graph 

We now demonstrate how the effect of hft-off described in Section 4.3 can create memory in 
a network. We use the decision graph with the simplex realisation as an example, and induce 
lift-off in the 2:3 direction at ^5. We give two examples: in the first the lift-off is maintained 
until the trajectory reaches ^2- In this example, we expect that trajectories which visit ^5 
will then have a larger probability of visiting ^3 on the next circuit of the network than 
they would have done had they previously visited ^q, ^7 or ^g- In the second example, we 
exhibit parameters such that the lift-off is not maintained at ^2, and consequently there are 
no memory effects. 

The particular structure of the decision graph can be used to simplify our analysis. Note 
that the sequence of equilibria visited (that is, the sequence ik, k = 1, . . .) can be deduced by 
only recording which of the equilibria ^5, ^6) and was visited by the trajectory on each 
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Figure 12: For the decision graph and simplex reahsation, the figure shows number of visits 
the trajectory makes to the equihbria ^4 (red), (blue) and (green), as a proportion of 
the number of visits to ^2, ^3 and ^4 respectively, for various noise levels (; the bars show 
standard errors. For details and parameter values, see text. 



'loop' around the network, and we therefore focus just on transitions between this subset of 
equilibria. 

Specifically, we consider the sub-itinerary l^. of (10) that consist of visits to {5, 6, 7, 8} by 
choosing the minimal strictly increasing set of indices < ik such that 

Ik = e {5,6,7,8}. 

Observe that from the structure of the decision graph in Figure 5 we expect that tk+i — ^k = 4 
on most occasions, and note also that we can have one-cycles (non-trivial transitions from j 
straight back to j) in this induced graph. Using the corresponding definitions for proportions 
of visits and transition probabilities between this subset of equilibria we define for any 
sequence of jk € {5, 6, 7, 8}: 

'P(ji) • • • ! Jm) = Prob( (a;, omega) : li = ji for i = 1, . . . , m, for the trajectory starting at {x, u)). 

(19) 

Analogously to (13) we can compute this from a typical trajectory as 

'P{ji, ■■■,jm) = lim 7#{0 < i < k : li+n = jn for = 1, . . . , m} (20) 
and so define the probability vfj of the first epoch being close to j G {5, 6, 7, 8}: 

= VU) 

and (assuming tij > 0) we define the transition probability between vertices and by: 

^iij2 = — ^ • (^1) 
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We perform numerical integrations using the Heun method for the cases with and with- 
out memory, and observe the effects in two ways. First, we compute the observed number of 
transitions between this subset of equihbria to estimate the matrix of probabihties 7ij-^j2 for 
ji,j2 = 5, ... ,8. A x-squared test can be used to determine whether or not the probabihties 
of visiting each of the equihbria are independent of the previous state of the trajectory. We 
also directly observe the distribution of coordinates as trajectories enter and leave neigh- 
bourhoods of equilibria. From the results in Section 4.3, we know that lift-off at ,^5 in the 2:3 
direction occurs if ^ < 1. This lift-off will still exist after the trajectory when the trajectory 
enters a neighbourhood of ^2 if ^ + < 1- We first perform an experiment where both 
conditions are satisfied, and then an experiment where only the first condition is satisfied. 



5.3 An example with memory 

We choose parameters so that 

^<1 and ^ + ^<1, 

esi 651 ei2 

and integrate the equations for the decision graph simplex realisation, with noise ( = 10~^, 
and for total time 200, 000, giving a sequence of 31, 568 epochs and so = 31, 568/4 cycles. 
The parameters used are Cjk = 2 for all j, k except for C13 = 0.8, C53 = 0.8, tjk = 2 for all j, 
k, ei2 = 2 , 624 = 636 = 643 = 1-9, 623 = 635 = 647 = 1-99, and 651 = cei = 671 = csi = 1.98. 

For a randomly chosen initial condition and noise path, the following matrix is the ob- 
served number of transitions between equilibria ^5, ^e, ^7 and in the reduced system: 



/2829 


1478 


143 


76 \ 


1166 


643 


379 


181 


355 


158 


88 


57 


\ 176 


90 


48 


25 J 



which gives an approximation to the matrix tt 



/ 7r5,5 


7r5,6 


7^5,7 


^5,8\ 




/0.6251 


0.3266 





0316 


0.0168\ 


7^6,5 


7^6,6 


7^6,7 


7r6,8 




0.4922 


0.2714 





1600 


0.0764 




7r7,6 


7r7,7 


7r7,8 




0.5395 


0.2401 





1337 


0.0866 


Vs,5 


7^8,6 


7^8,7 


^8,8/ 




1^0.5192 


0.2655 





1416 


0.0737/ 



where standard errors are of order 1/^/k^, i.e. ±0.01. Hence, if a trajectory visits ^5, it 
is more likely to go via ^3 (and hence to ,^5 or rather than ^7 or ^g) on its next circuit 
around the network. A ^-squared test confirms this by rejection the null hypothesis that the 
transition probabilities are independent of the previous equilibrium visited. 

Figure 13 demonstrates the hft-off by showing distributions of the xi and X3 coordinates 
as the trajectory passes from equilibria ^5 to ^1. We plot the coordinates as the trajectory 
crosses appropriate Poincare sections: H^, H^^^, and using h = 0.1. On H^, 

X3 is of order h, and the coordinate of interest is Xi. Figure 13(a) shows this to clearly be 
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(c) Distribution of on (d) Distribution of ^3 on iJ™* 

Figure 13: Distributions of coordinates for trajectories which pass from equihbrium ^3 to 
(^5 to ^1, in the decision graph. Parameters satisfy ff^ + < 1, and so hft-off in the 
coordinate is expected and can be seen in the non-Gaussian distributions. 

approximately Gaussian centered at zero. Figures 13 (b) to (d) show distributions of the 
coordinate on H^^^, if™ and H°^^. In each of these three cases, the distribution is clearly 
not Gaussian, but instead is the sum of two Gaussians shifted to the right and left of zero. 
This demonstrates the lift-off of the trajectory both in the positive and negative directions, 
on different loops around the network. 

The conditional distributions in Figure 14 clearly show the memory effect via a scatter 
plot of the X'i and 0:4 coordinates of the trajectory on conditional on the previously 

visited cycle. The conditional distribution peaks are clearly different as only those that went 
past ^6 at the last cycle display lift-off. 
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Figure 14: Scatter plot of 0:3 and X4 coordinates showing the condition distribution of trajec- 
tories as they exit from ^1 on for the same trajectory plotted in Figure 13. Trajectories 
which visited ^5 on the previous loop around the network are coloured red, while trajectories 
which previously visited ^e, ^7 and are coloured black, blue and green (note that the blue 
and green points are plotted under the black points, making them difficult to see). Lift-off in 
the X3 direction for those points that are coloured red can clearly be seen to cause memory. 



5.4 An example with no memory 

We do a second example with same parameters as in Section 5.3, except that C13 

•^53 ^ 1 1 , C53 Ci3 

— < 1 but 1 > 1. 

651 651 ei2 



1.3, so 



This means we expect to see lift-off in the X3 direction at ^5, but after the trajectory has 
passed ^1 this lift-off will be 'squashed' and hence the distribution of the X3 coordinate as the 
trajectory enters a neighbourhood of ^2 will be Gaussian again. Thus the network will have 
no memory. Again, we integrate for time 200, 000, which gives a sequence of kc = 30, 892/4 
cycles. 

For an example run, the observed number of transitions between equilibria .^5, ^q, ^7, and 
^8 is: 

/1686 957 660 348\ 
896 446 371 186 
697 328 254 142 
\ 373 167 136 76 / 
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Figure 15: Distributions of coordinates for trajectories crossing if™ and having pre- 

viously visited ^5. Parameters in this case satisfy — < 1 but — + — > 1. Lift-off in 

^ ^'J -J 651 651 ei2 

the X3 direction can be seen as the trajectory enters a neighbourhood of ^1 but it has been 
compressed by the time the trajectory leaves a neighbourhood of ^1. 

which gives the following approximation to tt^^ 



/^5,5 


7f5,6 


7^5,7 


^5,8\ 




/0.4618 


0.2621 


0.1808 


0.0953\ 


^6,5 


7f6,6 


7f6,7 


7r6,8 




0.4718 


0.2349 


0.1954 


0.0979 






7f7,7 


71"7,8 




0.4905 


0.2308 


0.1787 


0.0999 


\^8,5 




7f8,7 


^8,8/ 




^0.4960 


0.2221 


0.1809 


0.1011/ 



with standard errors of order ±0.01. In this case the x-squared test does not reject the null 
hypothesis of independence for this transition matrix. 

Figure 15 shows the distribution of the 0:3 coordinate (for those trajectories which pass 
through ^5) on if™ and ii°^*. It can clearly be seen that there is lift-off in the X3 direction 
before ^1, but on exiting, the distribution has returned to being approximately Gaussian 
with zero mean; this is shown as a scatter plot in Figure 16 where, in contrast to the case in 
Figure 14, the distributions do not appear to depend on previous epochs and so the system 
is memoryless. 



6 Discussion 

We have shown that, in principle, any finite directed graph can be realised as an embedded 
attracting robust heteroclinic network for a coupled cell system. In fact we do not even re- 
quire the network to be connected (respectively, for the components to be strongly connected) 
- in these cases the resulting network will be disconnected (respectively, the asymptotic dy- 
namics will be contained within a subnetwork). There will clearly be issues at vertices of 
high degree for outgoing edges - there must be nontrivial dynamics such as additional saddle 
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Figure 16: Scatter plot of 0:3 and X4 coordinates showing the condition distribution of trajec- 
tories as they exit from ^1 on for the same trajectory plotted in Figure 15. Trajectories 
which visited ^5 on the previous loop around the network are coloured red, while trajectories 
which previously visited ^q, ^7 and are coloured black, blue and green. The distributions 
of all four sets of trajectories appears to be identical, indicating that there is no memory. 



points involved in "separating" trajectories that go to different edges. We have not been 
able to prove all we would like to concerning the dynamics of the constructed networks. In 
particular we have not, to our satisfaction, been able to characterize the behaviour of com- 
plete chain recurrent set that contains the network for the simplex and cylinder realisations. 
Nonetheless, our numerical examples suggest that the dynamics is well enough understood 
to be useful as a method of designing dynamical systems. 

The networks will be robust to perturbations of the parameters up to the point where the 
network structure changes at a number of possible bifurcations. At one extreme we expect, 
suggested by [9], that the network will be destroyed (but replaced by an excitable network) 
if the saddles bifurcate to stable nodes, while at the other extreme a resonance bifurcation 
can cause the network to lose stability, resulting in bifurcation of approximately periodic or 
other attractors from the network. 

6.1 Residence times near vertices 

For a given typical trajectory we define Ik{j) to be the probability that the trajectory is near 
^j, that is: 

hU) = {le{l,...,k}:ti=j} and define iV,(j) = 
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to be the cardinality of Ik{j), i-e. the number of times a trajectory visits the equihbria ^j. 
This can be used to define the mean residence time of the trajectory close to the vertex 



In this paper we do not examine properties of tj except to remark that this is an interesting 
and nontrivial question [9]. It was shown in [37] that the mean passage time of a trajectory 
past an equilibrium (that is, the length of time during which the trajectory remains in a 
neighbourhood of an equilibrium) for a noisy heteroclinic cycle in the limit of low noise C ~^ 



where Xj is the largest expanding eigenvalue for equilibrium and h is the neighbourhood 
size. For a noisy heteroclinic network, the passage time will be similar, though if there are 
several expanding eigenvalues, it is possible that the transition time may depend on the exit 
route, and the mean residence times may encode the perturbations of the system to inputs 
with non-zero mean [40]. 

6.2 Relevance to computational systems 

The constructions presented here give flexible tools for designing coupled cell systems that 
realise finite-state computational systems and suggest new ways to adapt coupled cell systems 
systems to "learn" networks, by modifying parameters (cf [9]). An interesting question is 
whether either of the constructions can adapted to explain neural system computations that 
proceed on a dynamical basis. Simple network models of coupled neurons [6] can give rise 
to various structures of heteroclinic network in "phase space" . Can the construction here be 
used to improve this? 

In applications, it may also be desirable to have control over the statistical properties of 
how a trajectory moves around the heteroclinic network. We may also want to control the 
occurrence of memory effects in the network where transition probabilities depend on the 
recent history of the route taken around the network; the methods in Section 4 offers a route 
into doing this. 

There are many further questions that one could ask about the resulting designed net- 
works - these include, for example: To what extent can one design a network that has not 
only the given network structure, but also a specified set of average residence times and/or 
transition probabilities (with or without "memory" - namely dependence on the path up 
to the vertex)? This is likely to be an interesting and challenging problem where inclusion 
of anisotropic noise may important, as in [10]. Other questions concern the limits on the 
memory effects within such systems and the appearance of memory at larger noise levels. 

Finally, one aspect of the realisation methods here is that they are not very compact - 
the dimension of phase space scales linearly with the number of vertices (resp. edges) for 
the simplex (resp. cylinder) realisations. This could be a barrier to using these results as a 
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paradigm for neural computation where the encoding may be very dense. By contrast, the 
number of vertices in the "odd graph" networks that can be found in n = 2k + 1 globally 
coupled phase oscillators have n\/ {k\{k + l)\) vertices [8]. However, as previously highlighted, 
the latter networks are not easily usable for computation because they are very complex while 
also being strongly constrained by presence of a large number of symmetries. However, it 
would be interesting to find a way to robustly and stably realise a heteroclinic network in a 
"minimal" dimension network. 
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